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Abstract 

We examine stochastic temperature fluctuations of the cosmic background 
radiation (CBR) arising via the Sachs- Wolfe effect from gravitational wave 



O perturbations produced in the early universe. We consider spatially flat, per- 

q . turbed FRW models that begin with an inflationary phase, followed by a 

mixed phase containing both radiation and dust. The scale factor during the 
mixed phase takes the form a(rf) = cirj 2 + C2f] + C3, where q are constants. 
During the mixed phase the universe smoothly transforms from being radia- 



tion to dust dominated. We find analytic expressions for the graviton mode 
function during the mixed phase in terms of spheroidal wave functions. This 
mode function is used to find an analytic expression for the multipole mo- 
ments (af) of the two-point angular correlation function C(j) for the CBR 
anisotropy. The analytic expression for the multipole moments is written in 
terms of two integrals, which are evaluated numerically. The results are com- 
pared to multipoles calculated for models that are completely dust dominated 
at last-scattering. We find that the multipoles (af) of the CBR temperature 
perturbations for I > 10 are significantly larger for a universe that contains 
both radiation and dust at last-scattering. We compare our results with re- 
cent, similar numerical work and find good agreement. The spheroidal wave 
functions may have applications to other problems of cosmological interest. 

FACS numbers: 98.80.Cq, 98.80.C, 98.80.Es 
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I. INTRODUCTION 



This paper considers the effect of primordial gravitational waves on the cosmic back- 
ground radiation (CBR). We consider spatially flat, perturbed Friedman- Robertson- Walker 
(FRW) universes that begin with an early inflationary phase. As the universe rapidly ex- 
pands, perturbations of the spatial geometry that are local in origin (eg. thermal fluctu- 
ations) are quickly redshifted both in amplitude and wavelength. After sufficient inflation 
these perturbations are no longer visible to an observer; the only perturbations that remain 
visible within a Hubble sphere are quantum-mechanical zero-point fluctuations. Because 
these perturbations extend to arbitrarily high frequencies they can not be redshifted away. 
Since the only significant perturbations remaining after inflation are zero-point fluctuations, 
we assume that the initial state of the universe was the vacuum state appropriate to de Sit- 
ter space, containing only the quantum fluctuations and no additional excitations. As the 
universe continues to expand after inflation, these quantum fluctuations are redshifted to 
longer wavelengths and amplified; one may think of this in terms of particle (graviton) pro- 
duction (as we do), non-adiabatic amplification, or super-radiant scattering. In the present 
epoch the occupation numbers of the graviton modes are large. This is not surprising. It 
has been shown that one may treat the collective effects of primordial gravitational waves 
as being due to the presence of a stochastic background of classical gravitational waves, and 
since gravitons are bosons, such an interpretation is only possible if the occupation numbers 
are large. 

Sachs and Wolfe |l[ showed how gravitational wave perturbations result in CBR temper- 
ature anisotropy. As photons from the CBR propagate, the paths they follow are perturbed 
by the metric perturbation hij, which in this discussion is due entirely to primordial grav- 
itational waves. The energies of the photons are perturbed, which results in temperature 
fluctuations from point to point on the celestial sphere. These temperature fluctuations 
are usually characterized by the two-point angular correlation function C{^f) defined on the 
celestial sphere. Here 7 is the angle between two points on the sphere. Most often the two- 
point angular correlation function is expanded in terms of Legendre polynomials, and the 
expansion coefficients or multipole moments are calculated. For a derivation of the angular 
correlation function for spatially flat cosmologies see the recent paper by Allen and Koranda 
0. Henceforth we will assume that the reader is familiar with this paper, which contains a 
detailed review of previous work on this problem, a comprehensive discussion of the physical 
motivation, and a detailed and self-contained "first-principles" calculation. 

Early work on this subject |||],|],|| assumed that the universe was completely dust dom- 
inated at last-scattering when the CBR decoupled, as did recent work by White which 
presented a concise derivation of the formula for the multipole moments due to tensor pertur- 
bations. Also recently, Grishchuk has adapted the terminology and techniques of quantum 
optics to the analysis ||. Grishchuk stresses the importance of the phase correlations be- 
tween the modes of the metric perturbations. A similar analysis by Allen and Koranda 
H using standard "quantum field theory in curved space" techniques found equivalent re- 
sults. We also showed that the now standard formula given by Abbott and Wise [||] and 
Starobinsky || for the Z'th multipole moment is a long wavelength approximation to an 
exact formula. Both the work by Grishchuk, and by Allen and Koranda assumed that the 
universe was completely dust dominated at last-scattering. 
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The first authors to consider the CBR anisotropy (within the framework of the Sachs- 
Wolfe effect) for a universe that was not completely dust dominated at last-scattering were 
Turner, White, and Lidsey |J. They used a "transfer function" to express the solutions of the 
wave equation for the gravitational wave amplitude in terms of the standard long wavelength 
formula given by Abbott and Wise. They found that the standard formula, which assumes 
that the universe is completely dust dominated at last-scattering, consistently underestimates 
the contribution of gravitational waves to the CBR anisotropy. More recent work by Ng 
and Speliotopoulos [W\ used numerical methods to integrate the wave equation and found 
similar results. Neither of these studies used analytic expressions for the gravitational wave 
amplitude (or equivalently, the graviton mode function). 

Other work has been done which does not directly use the Sachs- Wolfe formula to calcu- 
late the anisotropy due to gravitational waves. Crittenden et. al. [[II] numerically evolved 
the photon distribution function using first-order perturbation theory of the general rela- 
tivistic Boltzmann equation for radiative transfer, and included a Thomson scattering source 
term. Dodelson, Knox and Kolb |12| have done a similar numerical analysis. Both found 
that the standard formula of Abbott and Wise is only accurate for small I multipole mo- 
ments, and consistently underestimates the contribution of the gravitational waves to the 
CBR anisotropy for higher I moments. 

In this paper we give the first correct analytic expression for the graviton mode function 
in a cosmology that transforms smoothly from being radiation to dust dominated. (We 
correct a minor error in earlier work by Sahni Jl3| and Nariai ]14|] which claimed to find an 
analytic expression for the mode function.) We use this analytic expression and the Sachs- 
Wolfe formula to find an analytic expression for the multipole moments (af) of the angular 
correlation function Ci^y). The analytic expression for the multipole moments is written 
in terms of two integrals; we use numerical methods to evaluate these integrals and report 
numerical values for the multipole moments. We compare our results with those mentioned 
above, and find good agreement. 

The paper is organized as follows. In section || we reproduce general expressions for 
the angular correlation function derived in ||, and explain how one uses these formulae to 
calculate the multipole moments for any spatially flat, inflationary cosmological model. In 
section |TXT| we introduce our cosmological model, which "begins" with an infinite de Sitter 
phase followed by a mixed phase that contains both radiation and dust. Early in the mixed 
phase the universe is radiation dominated; later it transforms smoothly from being radiation 
dominated to dust dominated. In section [TV] we solve the massless Klein-Gordon equation 
(or wave equation) for the graviton mode function. During the mixed phase the solutions 
to the wave equation are expressed in terms of spheroidal wave functions. The multipole 
moments are calculated in section [V| using the graviton mode function determined in section 



IV . An analytic expression for the moments is given in terms of two integrals, which are 



then evaluated numerically. In appendix A we discuss the spheroidal wave functions. The 
differential equation of spheroidal wave functions is introduced and its solutions examined. 
We introduce a useful notation for the spheroidal wave functions, and give a practical method 
for evaluating them. Finally in appendix B we describe the numerical techniques used to 
evaluate the spheroidal wave functions and the two multipole moment integrals. 

Throughout this paper we use units where the speed of light c = 1. We retain Newton's 
gravitational constant G and Planck's constant h explicitly. 
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II. THE ANGULAR CORRELATION FUNCTION 



We consider only the anisotropy of the cosmic background radiation (CBR) arising via 
the Sachs- Wolfe effect from tensor perturbations (gravitational waves). The anisotropy 
is characterized by the two-point angular correlation function C(7), where 7 is the angle be- 
tween two points located on the celestial sphere. The correlation function may be expanded 
in terms of Legendre Polynomials as 

C(7) = (f (0)^(7)) = E 2^(a?)Pz(c 0S7 ). (2.1) 

The expansion coefficients (af) are referred to as the multipole moments. The multipole 
moments are given in terms of an integral over graviton wavenumber k, 

(a, 2 ) = 4tt 2 (/ + 2){l + 1)1(1 - 1) jf T \k{k)\\ (2-2) 

where the function h(k) is proportional to the Sachs- Wolfe integral along null geodesies, 
and is given by 

w.pa^n faivl . (2.3) 

Jo k{r] - r] is - A) 2 

Here r/i s is the time of last-scattering, and rjo is the conformal time today. The function ji (z) 
is a spherical Bessel function of the first kind JTJJ. The function F(X, k) is proportional to 
the first derivative of the graviton mode function 0(r/, k), and is defined by 



F(X, k) = k 1/2 



^<P(v,k) 



»7=J7ls+^ 

The graviton mode function obeys the massless Klein-Gordon or wave equation 



(2.4) 



+ 24^0 + ^0 = 0, 2.5 



where a{rf) is the scale factor and 

drj 

The mode function must satisfy the Wronskian normalization condition 



(2.6) 



k)<j)*(r}, k) - P(rj, k)<j)(r], k)} = ~^- y (2.7) 

The only physical input required is the choice of an initial quantum state for the gravitational 
field, which amounts to a choice of boundary conditions for the wave equation ( |2.5| ) . 

These formula for the angular correlation function are very general; a detailed and com- 
plete derivation is given in [0. To calculate the correlation function (or equivalently the 
multipole moments) for any particular cosmological model, one need only to solve the wave 
equation (|2~5| ) for the graviton mode function, and substitute into the formulae ( |2.1| - [2.4D . 
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III. THE COSMOLOGICAL MODEL 

The spacetime considered here is a spatially flat, perturbed FRW universe. The metric 

ds 2 = a 2 (ri)[-dr] 2 + (<% + ^(77, x fc ))cfoW], (3.1) 

where 8ij is the flat metric of R 3 , rj is the conformal time, and a{rj) is the cosmological length 
scale or scale factor. The scale factor satisfies the Einstein equations 

where p(rj) is the energy density and P{rj) the pressure of the cosmological fluid. The metric 
perturbation hy is assumed to be small; in the limit as vanishes the spacetime is an 
unperturbed FRW universe. We have chosen a gauge so that the tensor perturbation 
has only spatial components. 

In the absence of the metric perturbation hij, the background spacetime is a spatially 
flat FRW universe. Since the spacetime is spatially flat, the density parameter Q , which 
is the ratio of the present-day energy density p to the critical energy density required to 
produce a spatially flat universe, is equal to unity: 

n ° = w = L (3 ' 3) 

To specify the model, we need to give the scale factor a(r]). We do this in such a way that 
the model is completely defined by the minimal set of free parameters given in Table |. All 
our results, including the final expression ( |5.2j ), can be expressed in terms of this minimal set 
of parameters. For clarity, we often define auxiliary quantities and express results in terms 
of them; the auxiliary quantities can be expressed in terms of the parameters in Table |. In 
typical inflationary models, the free parameters in Table | have values of order H between 
50 and 100 km s" 1 Mpc" 1 , 100 < Z h < 1500, 2 x 10 3 < Z cq < 2 x 10 4 , and 10 20 < Z end . For 



a review of inflationary cosmology see reference ]I6| 



A. The Inflationary or de Sitter Phase 

Our cosmological model passes through two phases. The first phase is a de Sitter or 
inflationary phase. During the de Sitter phase the universe expands exponentially (expressed 
in terms of comoving time t, the scale factor behaves like a(t) ~ e Ht , where H is the Hubble 
constant). In terms of conformal time, the scale factor is 

0(77) = a(r/end) (2 — ) for 77 < r/ end , (3.4) 

V ??end / 

where r/ en d is the conformal time at the end of the de Sitter phase. In terms of the parameters 
listed in Table |, 
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a(^end) = (1 + Z CQd ) 1 , 



(3.5) 



Vend 



[(2 + Z cq + Z end )(l + Z end )_ 



J cnd 



(3.6) 



where we have set the scale factor today a(i] ) = 1. During the de Sitter phase the energy 
density pas is constant and given by 



PdS 



3 d 2 (r] cnd ) 3H$(l + Z end ] 



87rGa 4 ( Vcnd ) 8vrG (2 + Z eq ) 



■(2 + Z cq + Z en d) 



- i7 2 ^ nd 
87rG 



'eq 



^end 



'eq 



Po- (3.7) 



The pressure during the de Sitter phase is negative and constant, and is given by 

-PdS = "PdS- 

The de Sitter phase is followed by a mixed phase. 



(3.8) 



B. The Mixed Radiation and Dust Phase 

Immediately following the de Sitter phase is a mixed phase containing both dust and 
radiation. The scale factor is 



a(r]) = a(r/e n d) 



7/ 



-1 + 



7] 



A VI + £/ \r] cnd / 7] end 
The constant £ is defined in terms of the free parameters by 

l + Z, 



for 1] > r] end . 



'eq 



1 + z end 

The stress-energy tensor is that of a perfect fluid with energy density 

P(r]) = Pdust(^) + Prad(^) for 7] > 7] cnd , 

where the energy density for the dust Pdust(^), and for the radiation p Tad (r]), are 

Pc q a 3 (r7 cq ) 



Pdust(^) = 
Prad(^) 



2 a 3 (?,) ' 

Peg a A {y cq ) 
2 a*( V ) ■ 



Here p cq = p(r) eq ) is the energy density at conformal time r] cq , when the dust and 
energy densities are equal, and is given by 



Peq 



3Hq (1 + Z eq j 1 ^ f ,, 



AnG (2 + Z, 



'end 



Pds ~ 2Z 3 p . 



(3.9) 

(3.10) 
(3.11) 

(3.12) 
(3.13) 
radiation 

(3.14) 



The pressure of the cosmological fluid P(rf) is 
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P(V) = = Pe< 



6q a\ V ) 

where the pressure at dust /radiation equality P eq is 



for 77 > 77 end , 



cq 



(3.15) 



(3.16) 



By inspection of (|3 . 1 3|) , (|5TT5|) , and fl3.16|) one sees that 



P(v) = Przd(ri) = -Ptadiv), 



(3.17) 



which is the relation between the pressure and energy density one would expect since the 
dust has zero pressure and the pressure is due entirely to the radiation present. The scale 
factor at r/ eq is 



cqj 



[1 + Z, 



cq ) 



(3.18) 



One can express r/ eq in terms of the conformal time at the end of the de Sitter phase 77 en d 
and the constant £ as 



Veq. 



Vend 



2^2(1 + 0-2-0^ 



^2(^2-1) 



(3.19) 



For 7] 7] eq the energy density (|3.11| ) varies as p(rj) ~ a _4 (^) and the cosmological model 
is radiation dominated, and for r\ 3> i] eq the energy density varies as p(r]) ~ a -3 (77) and 
the model is dust dominated. The model makes a smooth transition at rj = rj eq from being 
radiation to dust dominated. 



IV. THE GRAVITON MODE FUNCTION 

To calculate the multipole moments for the cosmological model described above, one 
must first solve the wave equation (|2.5| ) for the graviton mode function. We first solve for 
the "natural" positive- and negative-frequency mode functions during both the de Sitter and 
the mixed phases. We then make an appropriate choice of initial mode function during the de 
Sitter phase. This choice of initial mode function completely determines the graviton mode 
function for all later times. We express the mode function at later times using Bogolubov 
coefficient notation. 



A. The de Sitter Phase 



The solution to the wave equation fl2.5|) for the de Sitter phase can be expressed in terms 
of spherical Hankel functions. The scale factor during the de Sitter phase is given in ( |3.4| ). 
By making the change of dependent variable 

X (v,k) = (r ] -2r ]cnd )- 2 (j)(r ] ,k), (4.1) 



7 



and the change of independent variable 



z = k(r]-2r] cnd ), (4.2) 
the wave equation can be expressed in the form 

d? X . 2d X . I, 2 



«fa» + iS + l 1 -Fj x = a (4 ' 3) 

This is Bessel's differential equation, and the solutions are spherical Bessel or Hankel func- 
tions [15|. Using the normalization condition ( 2/7 ) one obtains for the positive-frequency 



mode function during the de Sitter phase 



where h'i\z) is a spherical Hankel function of the second kind |15|| , and pp = 1/hG 2 
5 x 10 93 gm/cm 3 is the Planck energy density The negative- frequency mode function dur- 
ing the de Sitter phase 4>l^\ri } k) is just the complex conjugate of the positive-frequency 
mode function ( |4.4|) . The positive- and negative- frequency mode functions form a complete 
solution to the wave equation for the de Sitter phase. 



B. The Mixed Radiation and Matter Phase 



The wave equation during the mixed radiation and dust phase (r/ > r/ C nd) can be cast 
in the form of the spheroidal wave function differential equation. By making the change of 
dependent variable 



X(v,k) = a 1/2 (r])(j)(r], k), 
the wave equation (|2~5| ) becomes 

_i_ v _i_ fu2 ld2 fa) ld fa) 



a{rf) \ Aa 2 (rj) 2 a(r)) 
One may define a new independent variable 



x = o. 



x 



Y 1 



i + 



so that 



dx 



1 a(rf) 



Using (|4.7| ) and 



<k] 2x a(r) cq ) 
one may write the wave equation (|4.6| ) in the form 



(4.5) 



(4.6) 



(4.7) 



(4.8) 
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fx 

dx 2 



+ 



2x 



a 2 (r]) 



1 

x 



2x 



Q(^cq) 

a{r]) 



+ 



4k x 



2J1 a2 (Veq) 



d 2 (r]) 



X 



2 Q 2 faq) 

a 2 {rj) 



2x 



dx 
dx 

2 d{rj)a 2 (r] C(l ) 



d 2 {ri)a{r]) 



X 



This expression can be simplified using the Einstein equations 
( |4.7|) One can show that 



d (77) = 2x 0(77)0(77, 



- 0. (4.9) 
along with ( |3.11| ) and 

(4.10) 



and using this expression one may write the wave equation 



as 



dx 2 

Since from ( |4. 7] ) 



d 2 Y d(Vea) d>Y 



0(77) dx 



2fc 



a(f], 



a(77 eq ) 2 a2 ( r 7eq) a(?7 eq ) 



x 



a 2 {rj) a(rj) 



X 



0. 



(4.11) 



a(rj) 



x 2 — 1 ' 



(4.12) 



and from Q and (^JJ 



0(?7eq) 



2^nd 



(4.13) 



the wave equation becomes 
d 2 x 2x 



dx 



dx 2 (1 — x 2 ) dx 
Here k is defined by 



+ 



(l-x 2 



+ 4k 



l-x 



2^2 



X 



0. 



(4.14) 



s c2 



t 2 2 
K ^end 



4tt 2 (2 + Z, 



eq y 



4tt z 



A^ 2 (1 + Z eq ) 2 



(4.15) 



eq 



where Ao is the present day wavelength of the mode with wavenumber k. For the multipole 
moments of interest, I G [2,1000], the contributions typically come from wavelengths in 
the range k G [10~ 3 , 10 3 ]. This form of the wave equation is the spheroidal wave function 
differential equation (|A1| ), which is discussed in detail in appendix A. 

The solutions to the wave equation ( |2.5| ) during the mixed phase can be expressed in 
terms of spheroidal wave functions. By inspection of (|4.7p one sees that x > 1, and so the 



possible solutions to the wave equation (|4.14j ) may be expressed as sums of any pair of 

X {x,k) = 3 SI(x,k), j = 1,2,3,4, (4.16) 



where J Sn( spheroidal wave function. Using (|4.5| ), ( [4.16|) , the normalization condi- 

tion ( p.7| ), and the Wronskian relation (|A22| ), one obtains a positive-frequency mode function 
during the mixed phase 
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&(V, k) = -2i, 



PdS a{n 



'end; 



\ 3vr p P a(r]) \ £ 



(4.17) 



Note that x is a function of the conformal time rj and k depends on the wavenumber k. The 



"*mixl 



(77, k) is just the complex conjugate of 



negative-frequency mode function ^[Jjj], k) 
the positive-frequency mode function. The positive and negative frequency mode functions 
form a complete solution to the wave equation during the mixed radiation and matter phase. 
Note that the choice of graviton mode function during the de Sitter phase completely deter- 
mines the mode function at all later times. Thus, the choice ( 4.17|) of "positive-frequency" 
during the mixed phase is unimportant. 



C. The Graviton Mode Function Expressed Using Bogolubov Coefficient Notation 



The choice of mode function during the de Sitter phase completely determines the mode 
function at all later times. This is because a solution to the wave equation depends only on 
the values of and on a spacelike Cauchy surface (ie. a surface of constant rf). We choose 
the mode function during the de Sitter phase to be the pure positive-frequency de Sitter 
solution C |4.4|) . This is the unique solution corresponding to a de Sitter-invariant vacuum 
state with the same (Hadamard) short distance behavior as one would find in Minkowski 
space [|17|- Having made this choice for the mode function during the de Sitter phase, the 
mode function during all times is 



(f)(rj,k) 



V < "^end, de Sitter phase, 



Q^end^inb^ 7 ?; &) + P (kVend) fimiLiV , k ) V > Vend, mixed phase 



(-), 



(4.18) 



where a and (3 are Bogolubov coefficients. 

The Bogolubov coefficients are determined by requiring that the mode function and its 
first derivative be continuous at 77 = 77 en d- One obtains 



a(kr) end ) 
P(ki] cnd ) 



kf/end 1 



'1+1 



vt+T 



A Sl'(xi, k) — 4 Sl(xi, k) 



vt+T 



2fc?7end 

i 



2kr], 



end 



1 + 



1 + 



kf]end 

i 



kr), 



end 



The prime in ( 4.19|) is defined by 
where 



(4.19) 
(4.20) 



Z=X\ 



a(r], 



end / 



1 + e 



(4.21) 
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Using the Wronskian relation for the spheroidal wave functions (|A22 ) one can easily verify 
that 



\a(kr) erid )\ 2 - \/3(kr] end )\ 2 = 1. (4-22) 

We stress again that the choice of graviton mode function during the de Sitter phase com- 
pletely determines the mode function at all later times. Thus, the choice of "positive- 
frequency" during the mixed phase is unimportant. Had we chosen a different solution 
to the wave equation during the mixed phase to call "positive-frequency" , then the Bogol- 
ubov coefficients ( [4.19| ) would be different in such a way so that the mode function (|4.1 
would be the same. 



V. MULTIPOLE MOMENTS 

Having determined the graviton mode function for our cosmological model, we may 
calculate the angular correlation function C(7), or equivalently, the multipole moments 
{of). We simply substitute the graviton mode function (|4.18|) into the formulae ( |2.1| - |2.4|) . 



A. Analytical Results 



The graviton mode function (|4.18j) is exact; no approximations (ie. long wavelength) 
have been made. One may directly substitute the mode function into the formulae (|2.1| - [2.4|) 
to obtain an exact expression for the multipole moments; because the arguments of the 
spheroidal wave functions in ( [4.18 ) are functions themselves, however, the result is compli- 
cated and not very illuminating. In typical inflationary models, the amount of expansion 
is very large. If one takes the limit Z end — > oo, then £ — > 0. This allows one to write a 
fairly compact expression for the multipole moments; after substituting the mode function 
Q4.18] ) into the formulae ( j2.1H2.4| ), one may collect together terms in the expression for the 
multipole moments (af) which are the same order in £. Then in the limit and 
£ — ► 0, one may consider only the leading term. The leading term, which is 0(£°), is given 
below in ( |5.2|) . The neglected terms are 0(£ x ) or greater. For typical inflationary models 
one has Z cnd > 10 20 and Z cq « 10 4 , so that £ < 10~ 16 . So the expression (|5.2p is quite 
accurate for most cosmological models, since the neglected terms are very small. 

The expressions for the multipole moments (af) are fairly simple. After substituting the 
mode function (|4.18| ) into the formulae (|2.1| - |2.4|) , and making the same changes of variable 



as in section II V H, ie. 



x 



\ X ^(^cq) 



and k 



*\ 1.2.. 2 



k r], 



end' 



(5.1) 



one obtains for the multipole moments (I > 2) 



(af) 



4 2 (/ + 2)!p dS 



-7T 



3 {l-2)\ppJo 



^{'TKxu «)G?(«) - 2 Tl( Xl , k)GI(k)} + 0(e)- (5.2) 
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The integral above is simply an integral over the (rescaled, dimensionless) wavenumber k. 
The functions ^T\(xi, k) are the same functions defined in (|A20| ). The functions G\{k) are 
the (reparametrized) Sachs- Wolfe integrals over null geodesies, given by 

px J; + i/ 2 (2/t 1 / 2 (x — x)) ( x . . v "I 

G\{k)= dx^ _ x)5/2 ^ 2 _ —^^—^S 2 (x, k) - >S 2 (x, «)|, (5.3) 

where the limits x\ s and x Q on the integral are given in terms of redshift by 



/ \ -\- Z 

x = J2 + Z eq and x is = Jl + eq . (5.4) 
v V 1 + Z ls 

These formulae, and especially the spheroidal wave functions, may be numerically evaluated 
using the techniques discussed in appendix B. Henceforth we refer to (|5.2|) and ( |5.3| ) as the 
"mixed" formulae for the multipole moments. Furthermore, we refer to multipole moments 
calculated using ( |5.2j ) and ( |5.3| ) as the "mixed" multipoles. 



B. Numerical Results 



We have numerically evaluated the "mixed" formulae (|5.2j ) and (|5.3| ). The second col- 
umn of Table [TI| shows the "mixed" multipoles. The third column of Table [H] shows the 
results obtained by evaluating Eq. (6.2) of Allen and Koranda ||, which we refer to as the 
"exact-dust" formula for the multipole moments. The "exact-dust" formula assumes (1) 
the universe begins with an initial de Sitter phase, followed by (2) a pure radiation phase 
containing only radiation and no dust, followed by (3) a pure dust phase containing only 
dust and no radiation, during which (4) last-scattering takes place. We refer to this type of 
universe as a "dust" universe, to distinguish it from the model described by (|3.9|), ( ft.llp , and 
( p.!5| ), which we refer to as "mixed" universe because it contains both dust and radiation 
at last-scattering. The fourth column shows the results obtained by evaluating Eq. (6.1) 
of Allen and Koranda, which we refer to as the "approximate-dust" formula for the multi- 
pole moments. The "approximate-dust" formula is a long-wavelength approximation to the 
"exact-dust" formula, and is equivalent to the standard formulae for the multipole moments 
given by Abbott and Wise f| and Starobinsky |J. The fifth column shows recent results 



of Ng and Speliotopoulos ||10|| , who numerically integrated the wave equation to obtain the 
amplitude of the gravitational wave (or graviton mode function). They also considered a 
universe model which contains both radiation and dust at last-scattering. All of the results 
in Table |I|were obtained with Z eq = 6000 and Z\ s = 1100. 

Figures [D,EL and 13] compare the "mixed" multipole moments to the multipole moments 



calculated using other techniques. The quantity M; in Figures |l]-0 is defined by 

(Note that in the equation defining M; contains an extraneous factor of Pds/pp-) The 
points labeled "mixed" in Figure [l] show results obtained by numerically evaluating the 
"mixed" formulae ( p.2[ ) and ( |5.3| ) . The points in Figure [I] labeled "transfer function" show 
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results from Turner, White, and Lidsey 0, who numerically integrate the wave equation, and 
express the solutions in terms of the standard long-wavelength approximate mode functions 
H using a "transfer function" . The points in Figure |2] labeled "mixed" again show results 
obtained by numerically evaluating the "mixed" formulae. The points labeled "exact-dust" 
show results obtained using the "exact-dust" formula, and the points labeled "approximate- 
dust" show results obtained using the "approximate-dust" formula. Figure |3| compares the 
"mixed" results to results from Dodelson, Knox, and Kolb [|TjJ (labeled "Boltzmann" ) , who 
do not use the Sachs- Wolfe formula to calculate the CBR anisotropy. Instead they use 
numerical methods to evolve the photon distribution function using first-order perturbation 
theory of the general relativistic Boltzmann equation for radiative transfer. All the multipole 
moments shown in Figures [l], and ||| are for cosmological parameters Z eq m 6000 and 
Z]a = 1100. (The values of Z eq for the Turner, White, and Lidsey work and for the the 
Dodelson, Knox, and Kolb results differs from 6000 by about one percent). 



C. Discussion 



The Ng and Speliotopoulos multipoles agree quite well with the "mixed" multipoles ob- 
tained here, using spheriodal wavef unctions, for a cosmology containing dust and radiation 
components. This is expected because the two methods used to calculate the multipole mo- 
ments should be essentially equivalent. The graviton mode functions Ng and Speliotopoulos 
obtain by numerically integrating the wave equation (with the correct boundary conditions) 
must be equivalent to our analytic expressions Q4.18|) . The cosmological model Ng and Spe- 
liotopoulos consider, however, is slightly different than our own. They model the smooth 
transition during the mixed phase by a simpler scale factor than our own ( |3.9| ). This may 
account for the small discrepancy between their results and our own. 

The Dodelson, Knox, and Kolb (Figure [|) multipoles also appear consistent with the 
"mixed" multipoles. This is also expected since the Boltzmann formalism |il~l , |l2l1 and the 
Sachs- Wolfe formalism should yield equivalent results. We do not understand at this time 
the small discrepancy between the two multipole spectrums. 

For I < 30, the "exact-dust" multipoles and the "approximate-dust" multipoles agree 
fairly well with the "mixed" multipoles (Figure [|). This is not surprising. The small I mul- 
tipoles (af ) are most affected by longer wavelength perturbations. These longer wavelength 
perturbations were redshifted outside the Hubble sphere early in the inflationary phase, and 
only recently re-entered the Hubble sphere (the longest wavelength perturbations remain 
outside the Hubble sphere even today [16|]). Because they remained outside the Hubble 
sphere until recently, the longer wavelength perturbations are insensitive to the details of 
the evolution of the universe before it became dust dominated. Thus, the "mixed" universe 
and "dust" universe models are essentially equivalent for longer wavelength perturbations. 
So it's not surprising that the same CBR anisotropies are produced by the long wavelength 
perturbations in either model. Thus the "exact- dust " , "approximate-dust", and "mixed" 
multipoles should, and do, agree for small I. Furthermore, because the "approximate-dust" 
formula for the multipole moments is equivalent to the standard formulae given by Abbott 
and Wise M and Starobinsky |J , one can conclude that the standard formulae for the mul- 
tipole moments are accurate for small I, whether or not the universe was completely dust 
dominated at last-scattering. 
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For I > 30, the "exact-dust" multipoles and the "approximate-dust" multipoles differ 
significantly from the "mixed" multipoles. Again this is not surprising. The larger I multi- 
poles are more affected by shorter wavelength perturbations, which re-entered the Hubble 
sphere before the universe became dust dominated, and are therefore sensitive to the details 
of the cosmological expansion before dust domination. A "mixed" universe becomes dust 
dominated much more slowly than a "dust" universe, which is dust dominated immediately 
after the radiation phase ends at dust/radiation equality r] e(l . Therefore the shorter wave- 
length modes in a "mixed" universe and a "dust" universe evolve very differently after they 
re-enter the Hubble sphere. This difference is evident in Figure |2] for the large I multipole 
moments. 

The Turner, White, and Lidsey (TWL) multipole spectrum differs significantly from the 
"mixed" multipole spectrum (Figure [I]). In particular their multipole moments are signif- 
icantly greater than the "mixed" multipoles for 3 < I < 80. Although TWL consider a 
"mixed" universe, they use the standard multipole moment formulae for a "dust" universe. 
However, they modify the standard formulae by including a time-independent transfer func- 
tion T(k/k eq ), which depends only on the wavenumber k. TWL give an explicit functional 
form for the transfer function, which they obtain by numerically integrating the wave equa- 
tion §: 

T(y) = [1.0 + 1.3% + 2.50y 2 ] 1/2 , where y = k/k eq . (5.6) 

Since the transfer function T > 1 appears in the expression for the multipole moments 
as \T{k/k eq )\ 2 (see Eqs. (22) and (23) of 0), one can see that the effect of the transfer 
function (|5.6j ) is to increase the multipole moments. Furthermore, the contribution from 
shorter wavelength (larger k) modes is enhanced more than the contribution from larger 
wavelength (smaller k) modes, since T(k/k eq ) — > 1 as k — ■> 0. 

Since the TWL multipoles are significantly greater than the "mixed" multipoles for 
3 ^ / ^ 80, the transfer junction overestimates the contribution to the multipole moments 
from longer wavelength modes. This is easy to see. As discussed above, the standard 
formulae are accurate for small I multipoles, and give essentially the same results as the 
"mixed" formulae for I < 30. The TWL formulae for the multipole moments are equivalent 
to the standard formulae, except for the transfer function. Thus if the transfer function 
was set to unity, the TWL multipoles would be the same as the standard multipoles, and 
hence would be equivalent to the "mixed" multipoles for I < 30. Since the TWL multipoles 
are larger than the "mixed" multipoles for 3 ^ / ^ 80, the transfer function must enhance 
the the contribution to the 3 ^ / ^ 80 multipole moments too much. Because the small / 
moments are affected most by longer wavelength perturbations, the transfer function must 
overestimate the contribution from longer wavelengths, or smaller wavenumber k. 

The TWL multipole spectrum is also significantly different than the "mixed" multipole 
spectrum for large I. The large I multipoles are most affected by shorter wavelength modes. 
The transfer function ( |5.6| ) significantly enhances the contribution of shorter wavelength 
modes (larger k) to the multipole moments. Because the TWL transfer function is time- 
independent, however, it can not alter the time evolution of the shorter wavelength modes. 
By comparing the "dust" universe graviton mode function (see Eq. (4.28) of 0) 



0Md „ l = ^___ 7 __ r (57) 
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(which is equivalent to the standard formulae for the gravity-wave amplitude @, and hence 
equivalent to the TWL formulae for the gravity-wave amplitude) to the "mixed" universe 
mode function (|4.18 ), one can see that the time evolution of the shorter wavelength modes 



is very different in a "dust" universe when compared to the time evolution in a "mixed" 
universe. Thus one expects that for shorter wavelength modes, the integral along the null 
geodesies ( [2.3|) in the Sachs- Wolfe formula will be very different for a "dust" universe, when 
compared to a "mixed" universe. Since the transfer function is time-independent, it has 
no effect on the integral (|2.3|) . So one expects that the multipole spectrum for a "dust" 
universe will be very different than the multipole spectrum for a "mixed" universe, even 
with the transfer function included. Figure [L] shows that this is the case. In particular, the 
"bump" in the TWL multipole spectrum appears near I rs 180, rather than I m 200, as for 
the "mixed" multipole spectrum. 



VI. CONCLUSION 

This paper examines the tensor perturbations of the gravitational field in a spatially flat, 
FRW cosmology containing a mixture of radiation and dust, and shows that they may be 
expressed in terms of spheroidal wave functions. Although spheroidal wave functions have 
appeared in this context before, previous authors incorrectly determined the characteristic 
exponent which labels these functions. After explaining the correct method for determining 
the characteristic exponent, we show that spheroidal wave functions may be efficiently and 
accurately evaluated using standard numerical techniques. 

We considered inflationary cosmological models, and used the spheroidal wave functions 
to find the spectrum of CBR temperature fluctuations resulting from primordial tensor 
(gravitational radiation) perturbations. These temperature fluctuations are predicted by 
all inflationary models. Their existence follows from first principles: it is a consequence 
of the uncertainty principle and the Einstein equation. The temperature fluctuations have 
been previously studied by a number of authors (including ourselves) using a variety of 
approximations, and both analytic and numerical techniques. 

In hindsight, only two approximations remain in this work. One is that the amplitude 
of the gravitational perturbation hij is very small. This approximation is indeed well jus- 
tified. The second approximation is that the energy density and pressure of the universe 
correspond to a mixture of dust and radiation as given in (|3.11| - [3T5l) . Going back in time, 
this approximation is good until approximately the time of nucleosynthesis, t = 200 seconds, 
when the number of effectively massless particles in the universe changed. This should not 
affect the multipole moments (of) which we consider, which have I < 300. For larger values 
of I, other physical effects such as the finite thickness of the surface of last scattering may 
also become relevant. 

It is likely that the spheroidal wave functions, which describe gravitational wave pertur- 
bations in realistic cosmological models, will find other useful applications. We expect that 
the results and methods of this paper may prove applicable to a wider variety of calculations 
than the CBR temperatures perturbations considered here. 
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APPENDIX A 

In spatially flat FRW universes, graviton mode functions in a linearized theory of gravity 



obey a minimally-coupled, massless, scalar wave equation |L8[ like (|2.5|) . If the scale factor of 
the spatially flat FRW universe transforms smoothly from being radiation to dust dominated, 
the solutions to the equation are spheroidal wave functions. 

1. The Differential Equation of Spheroidal Wave Functions and Its Solutions 

There is no generally accepted form for the differential equation of spheroidal wave 
functions. We write the differential equation as 

cP(p 2z dtp ( X An u 2 1 n 

" (T^jl + { (Trpj + 49 - JT^f /" = - (A1) 

The parameters /i, A, and 9 and the variable z can in general be complex. Here we take 
z,\,fx, and 9 to be real. We also consider only 9 > 0. The differential equation ( |A~1| ) has 
two regular singular points at z = ±1 and an irregular singular point at z = oo. We only 



consider z > 1. Then the solutions of (|A1|) are the spheroidal wave functions p9|j20 



<p = Srt\z,6), z>l, j = 1,2,3,4. (A2) 

The parameter [i is the order of the spheroidal wave function, and v is the characteristic 
exponent. In later sections we consider in detail the characteristic exponent v and its relation 
to the order \i and the parameters A and 9. For now we simply note that v is restricted so 
that 

u+ 1/2 ^ integer. (A3) 

For a very thorough and complete discussion of the differential equation ( |AI| ) and the solu- 
tions (|A2|) see reference . 

The spheroidal wave functions can be expressed in terms of more familiar special func- 
tions. If 9 = then the differential equation (|A1|) reduces to Legendre's differential equation, 



suggesting that the spheroidal wave functions can be expressed in terms of Legendre func- 



tions [19,20|. For z > 1, however, it is more useful to express the spheroidal wave functions 



as infinite sums of Bessel functions. For /i > one can write the spheroidal functions as 

-^-J T^\z,9), /i>0, (A4) 

where T^\z,9) is the infinite sum 
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T^\z,9) = s»(9) £ < r (0)vS 2 ,(20 1/2 2), J = 1,2,3,4. (A5) 

r=— oo 

The expansion coefficients a^ r (9) and the normalization factor s£(#), which are the same 
for any j, are discussed below. The functions ip^\z) are proportional to Bessel or Hankel 
functions: 



^>(*) = y/]fXil/2(z)- (A6) 

The spheroidal wave function of the first kind (j = 1) and the spheroidal wave function of 
the second kind (j = 2) form a complete solution to the differential equation ([All). Because 



the Hankel functions i^ 2 (2) = ± ii^(z) can be written as linear combinations of 



Bessel functions (see equation (3.86) in reference |15|]), the spheroidal functions of the third 
(j = 3) and fourth (j = 4) kind can be written as linear combinations of spheroidal functions 
of the first and second kind: 

s^\z,e) = s^\z,9)+ts^ 2 \z,e), 

S^(z,9) = S^(z,9) - iS^(z,9). (A7) 

The spheroidal functions of the third and fourth kind also form a complete solution to the 
spheroidal differential equation. The expansions ( |A4| ) of the spheroidal wave functions in 
terms of Bessel and Hankel functions are only useful if the infinite sums (|Al|) converge. 

The convergence of the infinite sums depends on the expansion coefficients a^ r {9). Sub- 
stituting ( |A4| ) and (|A5|) into the differential equation ( |A~1| ) yields a three-term recurrence 
relation that the expansion coefficients must satisfy. The recurrence relation can be written 
as 



where 



AU9) = 49 



A^ r (9)<r-i(0) + B£ r (*K >r (0) + C£ r (0)< r+1 (0) = 0, (A8) 
v + 2r - n)(u + 2r - ji - 1) 



(2z/ + 4r-3)(2z/ + 4r-l) 



BUM = A - {u + 2r)(u + 2r + 1) + (^+ 2r)(^ + 2r + 1) +^ - 1 
u ' rK ' V A ; (2z/ + 4r- l)(2z/ + 4r + 3) ' V ' 



(2z/ + 4r + 3)(2i/ + 4r + 5) 



The solution to this recurrence relation, as well as the convergence of the infinite sums 
(|A5|), depends critically on the parameters v, A, and 9, and is discussed below. Once a 
solution to the recurrence relation (for which the infinite sums (|A~5|) converge) is obtained, 
the normalization factor s%(9) is given by 
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-1 



(A10) 



This normalization is chosen so that in the limit as z becomes very large 



lim 



S^\z,9)/^\29 1/2 z) 



1. 



(All) 



This relation and many more details of the solutions (|A4|) are developed in reference |2(J 



2. The Eigenvalue A 



Although the order /i of the spheroidal wave function, along with the parameters 9 and 
A, appears directly in the differential equation (|A1|) , the characteristic exponent v does not. 
In most investigations and applications of spheroidal wave functions |19|,pO|,pT|,p2| , however, 



the parameter A is left unfixed; one assumes a (typically integer) value for v and considers 
A to be a function of /i, z/, and 9 and writes 



A = K{9)- 



(A12) 



\„{9) is often referred to as an eigenvalue, especially when considering spheroidal wave 



functions as solutions to the three-dimensional wave equation p3 |. 

For a given choice of the parameters /x, v, and 9, the eigenvalue \„{9) is that value of A 
for which the recurrence relation ( |A"8| ) has a minimal solution. A minimal solution, roughly 
speaking, is a set of coefficients a^ r {9) that satisfy the recurrence relation and fall off for large 
\r\ ||24j| . A dominant solution is a set of coefficients that satisfy the recurrence relation but do 
not fall off. If the solution to the recurrence relation is a dominant solution, the coefficients 



,(9) do not fall off for large |r|, and the infinite sums in (|A5[) may not converge. In 



general a three-term recurrence relation will have two independent solutions, much like a 
second-order, ordinary differential equation. However, neither of the two solutions, nor any 
linear combination of the two solutions, need be a minimal solution [24j|. For a given set of 



parameters /i, z/, and 9, a minimal solution to the recurrence relation ( |A8| ) exists only for a 
single, discrete value of A. That value for which the minimal solution exists is the eigenvalue 
\„{9). In our problem, we are given X,9, and fi. One can find v (modulo an integer) by 
requiring that the recurrence relation ( |A8| ) have a minimal solution. 

The functional relationship between the parameters n,v,9 and the eigenvalue \„{9) is 
complicated. It has been shown [|I^] that the functional relationship can be expressed as 

(A13) 



cos(27rz/) = /(A, [i A 



where a closed, analytic form for the function / is usually unattainable. The relation ( |A13| ) 
only determines the characteristic exponent v (as a function of A, fi, and 9) up to an integer; 
EDI fixes v: 



a second constraint 



0) = + 



(A14) 



In section [IV B| the special case of the differential equation ( |AI| ) with A = 2 is considered. 
The constraint ( |A14|) , along with the condition ( |A3| ) that v not be a half- integer, fixes v so 
that 
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1 3 

- < v < - tor A 

2 2 



2. 



(A15) 



For investigations of some of the analytic properties of /(A, /i 2 , (9) see reference [20|. A more 
practical method for finding the functional relation between the parameters /i, u, 9, and A is 
discussed in section |A~6. 



3. The Case fj, = 1, A = 2 



= 2 is of cosmological 
who incorrectly take 



The special case of the differential equation (Al) with // = 1 and A 
interest, and has been previously studied by Sahni [I3~| and Nariai [H 
the solution to be Sl^. In spatially flat FRW universes containing a mixture of dust and 
radiation, one may cast the wave equation for graviton mode functions in the form of the 
spheroidal wave function differential equation, as shown in section [IV Bj . The differential 
equation the graviton mode function obeys ( |4.14j) is the special case of the differential 

2, and 9 arbitrary. 



equation ( |A1[) with /i = 1, A = 2, and 9 arbitrary. [9 is arbitrary since in section [IV B| it 
is proportional to the wave number of a graviton mode function, and we desire expressions 
for the mode functions that are valid for an interesting range of wavenumbers.) Even if 
one takes A in to be fixed and does not consider the eigenvalue problem for A, the 
arguments and conclusion in the previous paragraph are still valid, especially the constraint 
( |A13| ); if fi, u, A, and 9 do not satisfy ( |A13| ) a minimal solution to the recurrence relation 
does not exist and the infinite sums in ( |A~5| ) do not converge. If A = 2, fi = 1, and 9 is 
arbitrary, then the characteristic exponent v is determined by the constraint ( |A13| ) with 
A = 2 and \i — 1. Tables of eigenvalues X^{9) for different values of fi, u, and 9 have been 
published [p2 |, however, and from these one can determine that the solution to ( |A13|) for 
/i = 1, A = 2, and 9 arbitrary is not v = 1, ie. 



cos(27r) 7^ /(2, 1, 9) for arbitrary 9. 



So for arbitrary 9 the solution to the differential equation 
the spheroidal wave function (|A2|) with /i = 1 and v — 1, ie 



for A = 2. 



(A16) 

l|) with /i = 1 and A = 2 is not 



(A17) 



This subtle point is missed in |Tj,[T4[, where the solution to the spheroidal wave function 



differential equation with /i = 1, A = 2, and 9 arbitrary is incorrectly given as S ± 



1(4) 



4. A More Suitable Notation 

The notation introduced above for the spheroidal wave functions (|A2|) is most useful 



when the order /i and the characteristic exponent v are fixed, and the eigenvalue A^(#) 
is considered as a function of /i, z/, and 9. Since we are most interested in the case when 
/i and A are fixed, it is more descriptive to denote the characteristic exponent as v = 
V\(9). This convention, however, would have one denote the spheroidal wave functions as 
S^\e){ z -i 9)i which is cluttered and inconvenient. In section |1V BJ where 9 is itself a function 
(of wavenumber), the notation would be even more unpleasant. 
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For this reason we adopt a new notation for the spheroidal wave functions. The solution 
to the differential equation (Al) with // and A fixed is denoted 



^ x (z,9) = S$ e) (z,9). (A18) 

The function on the left-hand side above is, of course, the same function as on the right-hand 
side, expressed using a different notation. The indices /i and A are those that appear in the 
differential equation, and again j = 1,2,3,4. In terms of Bessel functions, the spheroidal 
wave function is now written as 

^6)=(j^f' 2 n^6\ (A19) 

where 

oo 

n^z,e)=^{6) £ a^MMirW 1 ' 2 *), (A20) 

r=— oo 

with v = Vx{Q), and the ?/>^ are the same functions given in (|A6|) . The expansion coefficients 
a x (9) are the same as those in ( JA5|) and obey the same recurrence relation ( |A8| ) with the 
obvious change in notation. The normalization factor s x (9) is also the same as in (|A5|) with 



the obvious change in notation. We use this notation in sections |I V B[jVA| , and the rest of 
this appendix. 

5. Wronskian and Complex Conjugates 

Two relations for the spheroidal wave functions are especially useful in section [IV Bl The 



first involves the spheroidal wave functions of the third and fourth kind, and is obvious from 
the relations (|A7|) : 

*S£(M) = PSftM)]*- (A21) 

Here a * denotes complex conjugation. The second relation is the Wronskian for the 
spheroidal wave functions of the third and fourth kind: 

WfSftz, 9), 3 S^, 9)] ee 4 S^, 9) - 3 S^, 0)±*<%(z, 9) = gl/2( ^_ 1} - (A22) 



This relation, along with related results, is derived in reference |2tJ (see especially section 
3.65, equation (53) ). 



6. A Practical Method for Determining v x {6) and ^S x {z,6) 



The constraints ( [A13j ) and ( |A14| ) determine V\(9) for a given choice of //, A, and 9. In 
practice, however, it is difficult to use these constraints since a closed, analytic form for the 
the function / is usually not known. Recall, though, that the constraint ( |A13| ) is equivalent 
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to finding a minimal solution to the recurrence relation ( |A8|) . A minimal solution exists only 



for the single, discrete value of v^(Q) that satisfies the constraints ( A13 ) and ( A14Q . 

One can use the theory of continued fractions to find the minimal solution to the recur- 
rence relation and determine z^(#). Divide ( |A8|) by a r {9) and consider the infinite continued 
fraction Px r (9) defined for r > 1 as 



| 



U r+ i -+- U r +1" 



B r 



r+2 



where we have again suppressed for clarity the indices \x and A as well as the argument 9. 
For a given choice of parameters /i, A, 9, and v the infinite continued fraction can be used to 
find the ratios a r /a r _i of all the expansion coefficients, for r > 1, if the continued fraction 
converges. There is no loss of generality if one assumes ao = 1 |T9| , pO| , (note that this is 
compensated for in the normalization s^(9)) so the continued fraction (if it converges) can 
be used to find the expansion coefficients a r for r > 1. 

A second infinite continued fraction for r < — 1 can also be derived from the recurrence 
relation. We define the continued fraction 

KM - ±- = ~*V, = r < -1. (A24) 

Br + A ~ B ' + A ' -c r ; 



B r -2 + 



Assuming ao = 1, the continued fraction (if it converges) can be used to find the expansion 
coefficients a r for r < — 1. 

Whether or not the infinite continued fractions converge depends on the existence of 
minimal solutions to the recurrence relation. Pincherle's Theorem tells us that the 



infinite continued fraction P\(9) converges for r > 1 if and only if the recurrence relation 
has a minimal solution for r > 1. Further, if the infinite continued fraction does converge, 
it converges to a minimal solution. Likewise, N£(9) converges to a minimal solution for 
r < — 1 if and only if the recurrence relation has a minimal solution for r < — 1. 

A subtle, but important, point is that P\{9) and N£(9) may converge to different minimal 
solutions. Given a set of parameters /x, A, 9 and u, the continued fraction ( |A23| ) (with ao = 1) 

may converge and determine a set of coefficients ai, a^, These coefficients will satisfy 

the recurrence relation for r > 1 and will fall off for large r. Likewise, for the same set of 
parameters, the continued fraction ( [A 24] ) may converge and determine a set of coefficients 



a_i,a_2,a_ 3 , . . ., which will satisfy the recurrence relation for r < — 1 and will fall of for 
large \r\. The recurrence relation for r = 0, however, may not be satisfied. This is because 
the r = recurrence relation 

^,o(0K,-i(0) + KM<M + C%{9)a^{9) = (A25) 

is not explicitly solved when calculating either the P\{9) or the N^(9). One can see this by 
examining the continued fractions ( |A23| ) and (|A24|) ; the three coefficients a_i, ao, and a\ do 



not appear together anywhere in ( |A23| ) and ( |A24 ) for any r. So both the set of expansion 
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coefficients a%, a 2 , as, . . . found using P\{9) and the set <X—\ , 0_2, Ob— 3, • • • found using N%(9) 
are minimal solutions for a certain range of the index r, but not for all r. 

In order for the recurrence relation ( |A8| ) to be satisfied it must be true for all r. To 
ensure this, one must match the two solutions found using the continued fractions. This is 
accomplished by requiring that 

A%(9)N"x,-M + B^{0) + C%{6)P^{6) = Z»{9, v) = 0, (A26) 

which is just the requirement that the recurrence relation be satisfied for r = 0. So finding 
the minimal solution to the recurrence relation, and hence the characteristic exponent V\(9), 
is equivalent to finding the root of the the function Z%(9, v) defined in ( |A26| ). Given a test 
value for u^(9), one calculates the continued fractions P\i{9) and N£ and the rational 
functions A^ , B\ , and C^ , to find Z£(9, v). If Z% is not zero, one modifies the test value 
for Vx{9) using whichever root-finding algorithm one prefers. 

In practice, this method for determining v£(9) is efficient. Figure |] shows the function 
Z\{9, v) (of interest for the cosmological case) for 9 = 5 and ^ < v < §. This plot is typical 
for 9 in the range 10 -4 < 9 < 10 4 ; the root V\{9) is always located between two singularities, 
and Z\ is positive for u = u^{9) — e and negative for v = V\{9) + e where < e <C 1. This 
assists implementing a root-finding algorithm to find the zeroes of Z\[9, v) for arbitrary 9. 
Figure |5] shows the characteristic exponent v\{ff) for 10~ 3 < 9 < 10 3 . 

Once the characteristic exponent ^(#) is determined, one can use the continued fractions 
( |A23 ) and (|A24 ), along with a^ = 1 to calculate the remaining expansion coefficients. Since 



these coefficients are the minimal solution, the coefficients fall off for large |r|. One can then 
use (|A10 ) to find s^(9). Using ([A4Q and along with an algorithm for computing Bessel 



and Hankel functions, one can evaluate the spheroidal wave functions. Further details of 
the numerical techniques used may be found in Appendix B. Figure ^| shows the spheroidal 
wave functions ^(z, 9) and ^(z, 9) as functions of z for two values of 9. 



APPENDIX B 



This appendix briefly describes the numerical techniques used to obtain the results in 
section |VB| . One may separate the problem of numerically evaluating the spheroidal wave 
functions ^(z,^) and 2 S\{z,9) into two parts. The first is to determine the characteristic 
exponent v\(9) for a given 9. As noted in section |A 6| the characteristic exponent v\{9) is 



that value of v for which the function Z\{9, v) vanishes, where 

Z\{9, v) = A\ Q {9)Nl,-M + B\ fi {9) + C\ fi {9)PlM- (Bl) 

To find the root of Z^(9, v) (and hence the characteristic exponent v\{&) ) one must evaluate 
the five functions on the right-hand side of ( |B1| ). From (|A9|) one can see that the three 
functions A\ , B\ , and C\ (here A = 2, not v) are rational functions of v for any 9, and are 
easily evaluated. The continued fractions A 7 " 1 _ 1 (6') and Pii(9), defined by (|A23| ) and (|A24j ), 



are evaluated using the modified Lentz's method (see section 5.2 in reference [pq| ). Both 
continued fractions usually converge within 10 iterations when evaluated using this method. 
The root of Z\{9, v) is found using a simple bisection method. Although bisection may not 
be as efficient as other methods, it has the advantage that it is guaranteed to work once the 
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root has been bracketed. This is helpful since we are interested in finding the root v\{6) 
of Z\ (9, v) for many different 9; for some 9 the root v\{6) lies very close to a singularity of 
Z\{6, v), and in these cases other root-finding methods may not converge (see section 9.2 in 
reference f25|). 

Once the characteristic exponent v\(9) has been calculated, the remaining problem 
is to calculate the sums (|A20|) of expansion coefficients a\ r {9) times Bessel functions 
J v+2r+1/2 {29 l / 2 z) and Y v +2r+i/2(29 1 / 2 z). The expansion coefficients are calculated using 
the continued fractions Nl r {9) and P 2r {9). With a\ {9) = 1, one has 0^-^(9) = ^21^)1 
4,2^) = P2MPIM, and in general " 

n 

4nW = n P 2,W- (B2) 
J'=l 

Likewise, the negative index expansion coefficients are given by 

n 

4-»( fl ) = n jv 2,-iW- (B3) 

The continued fractions Nl r {9) and P 2r (9) are calculated using the modified Lentz's method 
as noted above. The normalization s\{9), defined in (|A10| ), is calculated at the same time 
as the expansion coefficients. 

The Bessel functions are calculated most efficiently using recurrence relations. The func- 
tions J v +2r+i/2{^9 l / 2 z) and Y u+2r+ i/ 2 {29 1 ! 2 z) with r > and r < are handled separately. 
For r > 0, J, +2r+1/2 (2^ 1 / 2 z) and J v+2r +\/2-2{^9 1 ' 2 z) are calculated for some large r using the 
routine bessjyQ found in reference P5| . Using the recurrence relations for Bessel functions 
(see 3.87 and 3.88 in |15[), the J u +2r+i/2(^9 1 ^ 2 z) are calculated using downward recursion to 
r = (this is the direction in which the recursion is stable for Bessel functions of the first 
kind |[25|| ). A similar procedure, using upward recursion, is used for the Bessel functions of 
the second kind Y v+2r+ i/ 2 {29 1 ^ 2 z). This gives the necessary Bessel functions for r > 0. 

The reflection formulae for Bessel functions are used to calculate the Bessel functions for 
r < 0. Since the index 

v + 2r + 1/2 = -(2|r| - v - 1/2) for r < and 1/2 < v < 3/2, (B4) 

the same procedure outlined above is used for J 2 \ r \~u-i/2('^9 1 ^ 2 z) and Y 2 \ r \_ u _i/ 2 {29 l l 2 z), and 
then the reflection formulae (see 6.7.19 of [ |25(| ) 

J-v(y) = cosvirJ v (y) - sin vnY v (y), (B5) 
Y-v(y) = sin vnJ v (y) + cos vnY v (y), (B6) 

are used to find Jv+ 2r +i/2(29 1 / 2 z) and Y v+2r+ i/ 2 (29 1 / 2 z) for r < 0. So all the Bessel functions 
needed are calculated with only a few time-consuming calls to the routine bessjyQ. 

Once the expansion coefficients a\ (0) and the Bessel functions J u +2r+i/2{^9 1 ^ 2 z) and 
^+2r+i/2(2^ 1 ^ 2 ^) are tabulated, the sums in ( |A20[ ) are calculated. Care must be taken when 



terminating the sums. Although the expansion coefficients a\ r {9) fall off very quickly as \r\ 
becomes large, the products a\ r {9) J u+2r+ i/ 2 (29 1 ^ 2 z) and a\^{9)Y y+2r+ i/ 2 {29 1 ^ 2 z) may not 
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fall off as fast. We terminate the sums when the products of the expansion coefficients and 
Bessel functions no longer contribute (at double-precision accuracy) to the sums. 

The primary numerical technique used to evaluate the multipole moments ( |5.2| ) is nu- 
merical integration. Both the integral over n in (|5.2| ) and the integral over x in (|5.3| ) were 
done using a fifth order embedded Runge-Kutta-Fehlberg algorithm with adaptive stepsize 
control ||25|| . Although formally the upper limit of the integral over k extends to infinity, 
we only integrated until the remaining contribution became negligible. The Bessel function 
with index I + 1/2 was evaluated with the routine bessjyQ. 
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FIGURES 



FIG. 1. Multipole moments (af) normalized to the quadrupole (a%). The horizontal axis is 
the index I of the multipole moment and the vertical axis is M;. See ( |5,5| ) for the definition of 
M/. The points labeled "mixed" show results obtained from the ( |5.2| ) and ( |5.3| ) which analytically 
model a universe containing both dust and radiation, and which is not completely dust dominated 
at last-scattering. The points labeled "transfer function" show results from Turner, White, and 
Lidsey || obtained using a transfer function. 



FIG. 2. Multipole moments (af) normalized to the quadrupole (a|)- The axes are the same 
as in Figure ||. The points labeled "mixed" show results obtained from (|5.2| ) and (|5.3| ) which 
analytically model a universe containing both dust and radiation, and which is not completely 
dust dominated at last-scattering. The points labeled "exact-dust" show results obtained using 
Eq. (6.2) of Allen and Koranda 0. The points labeled "approximate-dust" shows results obtained 
by evaluating Eq. (6.1) (which is a long-wavelength approximation to Eq. (6.2) ) of Allen and 
Koranda 0. Both sets of points "exact-dust" and "approximate-dust" are for universe models 
that are completely dust dominated at last-scattering. 



FIG. 3. Multipole moments (af) normalized to the quadrupole (ai,). The axes are the same 
as in Figure |l]. The points labeled "mixed" show results obtained from ( |5.2j ) and Q5.3| ) which 
analytically model a universe containing both dust and radiation, and which is not completely 
dust dominated at last-scattering. The points labeled "Boltzmann" show results from Dodelson, 
Knox, and Kolb |i"2fl , who do not use the Sachs- Wolfe formula to calculate the CBR anisotropy. 
Instead they use numerical methods to evolve the photon distribution function using first-order 
perturbation theory of the general relativistic Boltzmann equation for radiative transfer. 



FIG. 4. The function Z\{Q, v) plotted versus v for = 5. The value of v at which this function 
vanishes is the characteristic exponent v\(0) for = 5. This plot is typical for values of in the 
range 10~ 4 < < 10 4 . The root is always located between two singularities. 



FIG. 5. The characteristic exponent v\{Q) plotted vs log (9. Only for certain discrete values of 

(z, 0) is not a valid solution to the differential 



is v equal to 1. For this reason the function S 1 
equation (Al) for /i = 1, A = 2, and arbitrary 



FIG. 6. The spheroidal wave functions 1 S\(z,0) (solid lines) and 2 S\(z,0) (dotted lines) for 
= 0.1 and = 10. Both functions diverge in the limit as z — > 1. Note that the vertical scales are 
different for the two plots. 
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TABLES 



TABLE I. List of free parameters that define the cosmological model 



Parameter 



Units 



Range 



Description 



Z\ K 

Z C q 
Z C nd 



length -1 
dimensionless 
dimensionless 
dimensionless 



H >0 
Z ls >0 
unrestricted 

Zend > Z e q, Z\ s 



Present-day Hubble expansion rate 

Redshift at last scattering of CBR 

Redshift at equal matter /radiation energy density 

Redshift at end of de Sitter inflation 



TABLE II. Multipole moments (of) evaluated using different methods. These have been di- 
vided by the scale of the moments Pds/Pp- The second column shows multipoles for a mixed 
cosmology obtained from the formulae ( |5.2| ) and ( |5.3|) , which analytically model a universe con- 
taining both dust and radiation, and which is not completely dust dominated at last-scattering. 
The third column shows results obtained by evaluating Eq. (6.2) of Allen and Koranda Q, which 
assumes that the universe was completely dust dominated at last-scattering. The fourth column 
shows results obtained by evaluating Eq. (6.1) of Allen and Koranda, which is a long wavelength 
approximation to their Eq. (6.2), and equivalent to standard formula of Abbott and Wise [|| and 
Starobinsky 0. The fifth column shows results obtained by Ng and Speliotopoulos who nu- 
merically integrated the wave equation to obtain the amplitude of the gravitational waves. They 
also considered a universe model which is not completely dust dominated at last-scattering. All 
the results were obtained with Z cq = 6000 and Z\ s = 1100. 
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